% Having calibrated the firm side (incl. Hotelling parameter xi),
% this script does the following:
%
% 1) calculates actual profits for each Rz,
% 2) finds a Pareto-efficient contract (Rz that maximises U_0 s.t. staying
% on the same iso-profit),
% 3) calculates the associated CE welfare loss and the impact on savings
% outcomes by comparing the profit-maximising with the efficient contract,
% 03/2024 UPDATE:
% 3') calculates the Ramsey welfare loss (as if delta=1),
% 4) plots the consumption paths under different contracts. 


clear
load('recalibrated_data_with_eval.mat')

fileID = fopen('res_Firm_POSTcalibration_Analysis_2021.txt','wt')

rho = mydata.rho;
delta = mydata.delta;
effhh_size_ = mydata.effhh_size;
avgIncome = mydata.avgIncome;
meanhhs = mydata.meanhhs;
R = mydata.R;
zpen_eval = mydata.zpen_eval;
alpha = mydata.alpha;
survival = mydata.survival;
age_ = mydata.age;


THETAtarget = 0.05;
totalcosttarget = 0.005;
adminratiotarget = 0.50;
markuptarget = 0.2;

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_naif');
avgZ_naif_target = mydata.(suffix);
suffix = strcat('fmax_',num2str(THETAtarget*10000));
fmax_target = mydata.(suffix);


% Give me the calibrated values:
gamma1 = 0.5;
c1 = adminratiotarget*totalcosttarget*mean(avgZ_naif_target) / mean(avgZ_naif_target.^(gamma1));
gamma2 = 5.75;
c2 = (1-adminratiotarget)*totalcosttarget*mean(avgZ_naif_target) / (1+THETAtarget-R)^(gamma2);
phi = 0.35;

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['gamma1 = %f, c1 = %f, gamma2 =  %f, c2 = %f \r\n'],[gamma1, c1, gamma2, c2]);
fprintf(fileID, ['**************\r\n']);


% 1) Back out the Hotelling parameter xi
xifound = 0;

% 1a) Calculate utility from the (supposed) equilibrium contract
if rho == 1
    yearly_nobeq = meanhhs * log(avgIncome/meanhhs);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
else
    yearly_nobeq = meanhhs * ((avgIncome / meanhhs)^(1-rho)-1) / (1-rho);
    utility_nobeq = (yearly_nobeq / (1-delta)) - yearly_nobeq;
end

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_TC');
avgZ_TC = mydata.(suffix);
suffix = strcat('avgX',num2str(THETAtarget*10000),'_TC');
avgX_TC = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETAtarget*10000),'_TC');
avgTC_TC = mydata.(suffix);
    
bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_TC(t)+(1-zpen_eval(t))*avgZ_TC(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

avgTC_fee = avgTC_TC - phi*fmax_target*ones(1,71);
    
instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end
    
lifetime_equilibrium = instant(1);
for t = 2:71
    lifetime_equilibrium = lifetime_equilibrium + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

% 1b) Check that the firms don't want to deviate to lower/higher fees for target Rz
gridlength = 11;
gridjump = 0.05;

ProfitPVvec = zeros(1,gridlength);
expected_profits = zeros(1,gridlength);
xi = 0.0001;

for j = 1:gridlength
    
    % First, profit from an interaction
    fee = (1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*phi*fmax_target;
    
    THETA = THETAtarget;
    avgZ_ = avgZ_naif_target;
    
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);

    profits = ones(1,71)*fee - costs;

    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
    
    ProfitPVvec(j) = ProfitPV;
    
    % Second, probability of interaction
    avgTC_fee = avgTC_TC - fee*ones(1,71);
    
    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    lifetime_fee = instant(1);
    for t = 2:71
        lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
    
    probabilityraw = 0.5 + (lifetime_fee - lifetime_equilibrium)/(2*xi);
    probability = max(0,min(1,probabilityraw));
    
    expected_profits(j) = ProfitPV*probability;
end

if max(expected_profits) == expected_profits(ceil(gridlength/2))
    xifound = 1;
    fprintf(fileID, ['**************\r\n']);
    fprintf(fileID, ['For xi = %f, the calibrated fee of \r\n is Hotelling equilibrium'],[xi,phi*fmax_target]);
else
    xi = 1.25*xi;
end

counter = 1;
while xifound < 1 && counter < 1000
    expected_profits = zeros(1,gridlength);
    for j = 1:gridlength
        
        fee = (1 - floor(gridlength/2)*gridjump + (j-1)*gridjump)*phi*fmax_target;
    
        avgTC_fee = avgTC_TC - fee*ones(1,71);
    
        instant = zeros(1,71);
        for t = 1:71
            if rho == 1
                instant(t) = effhh_size_(t) * log(avgTC_fee(t)/effhh_size_(t));
            else
                instant(t) = effhh_size_(t) * ((avgTC_fee(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
            end
        end
    
        lifetime_fee = instant(1);
        for t = 2:71
            lifetime_fee = lifetime_fee + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
        end
    
        probabilityraw = 0.5 + (lifetime_fee - lifetime_equilibrium)/(2*xi);
        probability = max(0,min(1,probabilityraw));
    
        expected_profits(j) = ProfitPVvec(j)*probability;
    end

    if max(expected_profits) == expected_profits(ceil(gridlength/2))
        xifound = 1;
        fprintf(fileID, ['**************\r\n']);
        fprintf(fileID, ['For xi = %f, the calibrated fee of %f is Hotelling equilibrium \r\n'],[xi,phi*fmax_target]);
    else
        xi = 1.25*xi;
    end
    
    counter = counter + 1;
end


% 1c) Profit level at the target
THETA = THETAtarget;
avgZ_ = avgZ_naif_target;
    
costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);

profits = ones(1,71)*phi*fmax_target - costs;

maxprofit = profits(1);
for t = 2:71
    maxprofit = maxprofit + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
end


% 1d) Need to re-calculate the iso-fees
for THETA = mydata.THETA_grid
    if abs(THETA - THETAtarget) < eps
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = phi*fmax_target;
    else
        suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
        avgZ_ = mydata.(suffix);
                        
        costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
                        
        fee = 100;
        profits = ones(1,71)*fee - costs;
        ProfitPV = profits(1);
        for t = 2:71
            ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
        end
                    
        while maxprofit > ProfitPV
            fee = fee + 6;
            profits = ones(1,71)*fee - costs;
            ProfitPV = profits(1);
            for t = 2:71
                ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
            end
        end
                        
        suffix = strcat('fmax_',num2str(THETA*10000));
        fmax = mydata.(suffix);
                        
        suffix = strcat('isofee_',num2str(THETA*10000));
        workspace.(suffix) = min(fee-3,fmax);
    end
end

                    
% 1e) Calculate profits for each Rz (just to double-check the iso-profit condition)
fprintf(fileID, ['**************\r\n']);
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    
    costs = c1*avgZ_.^(gamma1) + ones(1,71)*c2*(1+THETA-R)^(gamma2);
    
    % This is just needed for calculating out-of-eqm predictions
    if abs(THETA - THETAtarget) < eps
        costs500 = costs;
    end
    
    suffix = strcat('isofee_',num2str(THETA*10000));
    fee = workspace.(suffix);
         
    profits = ones(1,71)*fee - costs;
    ProfitPV = profits(1);
    for t = 2:71
        ProfitPV = ProfitPV + (1/R)^(t-1) * prod(survival(1:t)) * profits(t);
    end
             
    fprintf(fileID, ['For Rz = %f, the PV of profits (on the iso-profit curve) is %f \r\n'],[THETA,ProfitPV]);
end


% 1f) Calculate which of the contracts lying on the iso-profit maximises
% ACTUAL consumer welfare. 
fprintf(fileID, ['**************\r\n']);

actual_lifetime_iso_grid = [];
for THETA = mydata.THETA_grid
    suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
    avgZ_ = mydata.(suffix);
    suffix = strcat('avgX',num2str(THETA*10000),'_naif');
    avgX_ = mydata.(suffix);
    suffix = strcat('avgTC',num2str(THETA*10000),'_naif');
    avgTC_ = mydata.(suffix);
                    
                    
    suffix = strcat('isofee_',num2str(THETA*10000));
    isofee = workspace.(suffix);
    avgTC_ = avgTC_ - ones(1,71)*isofee;

    instant = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant(t) = effhh_size_(t) * log(avgTC_(t)/effhh_size_(t));
        else
            instant(t) = effhh_size_(t) * ((avgTC_(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    bequest = zeros(1,71);
    for t = 1:71
        cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
        if rho == 1
            yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        else
            yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
            utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
        end
        bequest(t) = alpha*(utility_withbeq - utility_nobeq);
    end
    
    lifetime = instant(1);
    for t = 2:71
        lifetime = lifetime + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
    end
                    
    actual_lifetime_iso_grid = [actual_lifetime_iso_grid lifetime];
    fprintf(fileID, ['For Rz = %f and the associated iso-fee, the ACTUAL lifetime utility is %f \r\n'],[THETA,lifetime]);
                    
    if abs(THETA - THETAtarget) < eps
        actual_lifetime_at_target = lifetime;
    end
end


% Comparison between equilibrium & efficient contracts
[eff_benchmark, id_efficient] = max(actual_lifetime_iso_grid);

fprintf(fileID, ['**************\r\n']);
fprintf(fileID, ['The Pareto-efficient contract has Rz = %f \r\n'],[mydata.THETA_grid(id_efficient)]);


% What is the CE welfare loss?
utility_multiplied = actual_lifetime_at_target;

suffix = strcat('avgZ',num2str(THETAtarget*10000),'_naif');
avgZ_ = mydata.(suffix);
suffix = strcat('avgX',num2str(THETAtarget*10000),'_naif');
avgX_ = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETAtarget*10000),'_naif');
avgTC_ = mydata.(suffix);
suffix = strcat('isofee_',num2str(THETAtarget*10000));
isofee = workspace.(suffix);
    
avgTC_ = avgTC_ - ones(1,71)*isofee;

bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end

efficient = eff_benchmark;
multiplier = 1.00;
while utility_multiplied < efficient 
    instant_multiplied = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant_multiplied(t) = effhh_size_(t) * log((multiplier*avgTC_(t))/effhh_size_(t));
        else
            instant_multiplied(t) = effhh_size_(t) * (((multiplier*avgTC_(t))/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    utility_multiplied = instant_multiplied(1);
    for t = 2:71
        utility_multiplied = utility_multiplied + delta^(t-1) * prod(survival(1:t-1)) * (survival(t)*instant_multiplied(t) + (1-survival(t))*bequest(t));
    end
    multiplier = multiplier + 0.0001;
end

fprintf(fileID, ['The associated CE welfare loss is %f percent\r\n'],[100*(multiplier-1.0001)]);


% ***************************************************
% 03/2024: What is the Ramsey welfare loss (delta=1)?
% ***************************************************
THETA = mydata.THETA_grid(id_efficient);

suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
avgZ_ = mydata.(suffix);
suffix = strcat('avgX',num2str(THETA*10000),'_naif');
avgX_ = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETA*10000),'_naif');
avgTC_ = mydata.(suffix);
                                        
suffix = strcat('isofee_',num2str(THETA*10000));
isofee = workspace.(suffix);
avgTC_ = avgTC_ - ones(1,71)*isofee;

instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(avgTC_(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((avgTC_(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end
    
bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end
    
lifetime = instant(1);
for t = 2:71
    lifetime = lifetime + 1.0 * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

actual_lifetime_efficient_RAMSEY = lifetime;


THETA = THETAtarget;

suffix = strcat('avgZ',num2str(THETA*10000),'_naif');
avgZ_ = mydata.(suffix);
suffix = strcat('avgX',num2str(THETA*10000),'_naif');
avgX_ = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETA*10000),'_naif');
avgTC_ = mydata.(suffix);
                    
suffix = strcat('isofee_',num2str(THETA*10000));
isofee = workspace.(suffix);
avgTC_ = avgTC_ - ones(1,71)*isofee;

instant = zeros(1,71);
for t = 1:71
    if rho == 1
        instant(t) = effhh_size_(t) * log(avgTC_(t)/effhh_size_(t));
    else
        instant(t) = effhh_size_(t) * ((avgTC_(t)/effhh_size_(t))^(1-rho)-1) / (1-rho);
    end
end
    
bequest = zeros(1,71);
for t = 1:71
    cons_bequest = (R-1)*(max(0,avgX_(t)+(1-zpen_eval(t))*avgZ_(t))) + avgIncome;
    if rho == 1
        yearly_withbeq = meanhhs * log(cons_bequest/meanhhs);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    else
        yearly_withbeq = meanhhs * ((cons_bequest / meanhhs)^(1-rho)-1) / (1-rho);
        utility_withbeq = (yearly_withbeq / (1-delta)) - yearly_withbeq;
    end
    bequest(t) = alpha*(utility_withbeq - utility_nobeq);
end
    
lifetime = instant(1);
for t = 2:71
    lifetime = lifetime + 1.0 * prod(survival(1:t-1)) * (survival(t)*instant(t) + (1-survival(t))*bequest(t));
end

actual_lifetime_at_target_RAMSEY = lifetime;


utility_multiplied = actual_lifetime_at_target_RAMSEY;
efficient = actual_lifetime_efficient_RAMSEY;
multiplier = 1.00;
while utility_multiplied < efficient 
    instant_multiplied = zeros(1,71);
    for t = 1:71
        if rho == 1
            instant_multiplied(t) = effhh_size_(t) * log((multiplier*avgTC_(t))/effhh_size_(t));
        else
            instant_multiplied(t) = effhh_size_(t) * (((multiplier*avgTC_(t))/effhh_size_(t))^(1-rho)-1) / (1-rho);
        end
    end
    
    utility_multiplied = instant_multiplied(1);
    for t = 2:71
        utility_multiplied = utility_multiplied + 1.0 * prod(survival(1:t-1)) * (survival(t)*instant_multiplied(t) + (1-survival(t))*bequest(t));
    end
    multiplier = multiplier + 0.0001;
end

fprintf(fileID, ['(The associated CE welfare loss as measured by Ramsey is %f percent)\r\n'],[100*(multiplier-1.0001)]);
                    
% ***************************************************
% 03/2024: New code ends here.
% ***************************************************                    


% What is the impact on savings?
suffix = strcat('avgZ',num2str(THETAtarget*10000),'_naif');
avgZ_target = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETAtarget*10000),'_naif');
avgTC_target = mydata.(suffix);
suffix = strcat('isofee_',num2str(THETAtarget*10000));
isofee_target = workspace.(suffix);

THETAefficient = mydata.THETA_grid(id_efficient);
suffix = strcat('avgZ',num2str(THETAefficient*10000),'_naif');
avgZ_efficient = mydata.(suffix);
suffix = strcat('avgTC',num2str(THETAefficient*10000),'_naif');
avgTC_efficient = mydata.(suffix);
suffix = strcat('isofee_',num2str(THETAefficient*10000));
isofee_efficient = workspace.(suffix);

% t=45 corresponds to the age of 64 - final year prior to retirement
savings_impact = (avgZ_efficient(45) - avgZ_target(45))/ avgZ_target(45);
fprintf(fileID, ['Savings at retirement would be %f percent higher under an efficient contract \r\n'],[100*(savings_impact)]);


% Finally, plot the consumption paths elicited by optimal vs.
% efficient contracts (Figure 3). Multiply by 1.905 to express monetary
% amounts in 2018 USD, rather than 1990 USD.  

TCpath_target = (avgTC_target - isofee_target*ones(1,71))*1.905;
TCpath_efficient = (avgTC_efficient - isofee_efficient*ones(1,71))*1.905;

retirement_cons_impact = (mean(TCpath_efficient(45:end)) - mean(TCpath_target(45:end)))/ mean(TCpath_target(45:end));
fprintf(fileID, ['Consumption in retirement would be %f percent higher under an efficient contract \r\n'],[100*(retirement_cons_impact)]);


figure
plot(age_(1:end-10),TCpath_efficient(1:end-10),'b', 'LineWidth',1.0)
hold on
plot(age_(1:end-10),TCpath_target(1:end-10),'r--','LineWidth',1.25)
grid on
title(['Consumption path over the life cycle ($c_1 =$',num2str(c1, '%1.3f'),', $\gamma_1 =$ ', num2str(gamma1),', $c_2 = $', num2str(c2,4), ', $\gamma_2 = $ ', num2str(gamma2,4),')'],'interpreter','latex')
xlabel('Age')
ylabel('Total consumption (2018 USD)')
legend({['Total consumption under the Pareto-efficient contract ($R^Z =$ ',num2str(THETAefficient*100),'$\%$)'] ,['Total consumption under the profit-maximising contract ($R^Z = $ ',num2str(THETAtarget*100),'$\%$)']},'interpreter','latex','location','southeast')

   